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ABSTRACT 



By means of a one-zone evolutionary model, we study the co-evolution of supermassive black holes and their host galaxies, as a 
function of the accretion radiative efficiency, dark matter content, and cosmological infall of gas. In particular, the radiation feedback 
is computed by using the self-regulated Bondi accretion. The models are characterized by strong oscillations when the galaxy is 
in the AGN state with a high accretion luminosity. We found that these one-zone models are able to reproduce two important 
phases of galaxy evolution, namely an obscured-cold phase when the bulk of star formation and black hole accretion occur, and the 
following quiescent hot phase in which accretion remains highly sub-Eddington. A Compton-thick phase is also found in almost all 
models, associated with the cold phase. An exploration of the parameter space reveals that the closest agreement with the present-day 
Magorrian relation is obtained, independently of the dark matter halo mass, for galaxies with a low-mass seed black hole, and the 
accretion radiative efficiency ^0.1. 
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1. Introduction 

Elliptical galaxies invariably contain central supermassive black 
holes (SMBHs), and there exists a tight relationship between 
the characteristic stellar velocity dispersion cr (or stellar mass 
M„) of the host system and the SMBH mass M B h (e.g.. 



Magorrian et al.ll998t[Ferrarese & Merrittl200"ol: Tremaine et alj 
20021 Yu & Tremaine 2002). These relations clearly indicate a 



co-evolution of the SMBHs and their host spheroids. Several in- 
vestigations have been dedicated to th is subject, either by us- 
ing h y drody namical simulations (e.g.. ICiotti & Ostrikerl 1 1 997L 
2001. 120071: ICiotti et al I I2009L |201(J) or one-zone models 



(e.g., Sazonov et al. 2005 



.. _ hereafter SOCS;JMeio_etalj[2008, 
lMatteuccil2008l) . In some work, the effect of galaxy merging ha s 
also been taken into account (e.g., iHopkins et aTT l2005l 12006). 
The main advantage of hydrodynamical models is that complex 
physical phenomena effects (such as shock waves, jets, radiation 
transport, etc) can be taken into account. However, the computa- 
tional time of the simulations force us to search for faster meth- 
ods that allow a more systematic exploration of the parameter 
space, which is prohibitive with hydrodynamical simulations. In 
this framework, hydro-simulations can be used to set the accept- 
able range for parameters to be adopted in toy models (e.g., the 
duty cycle value, see Sect. 12. U . while simpler models, such as 
that used here, are useful to identify the most interesting cases 
that can be simulated in detail with hydrodynamical codes. 

The general idea behind one-zone models is to work with 
"average" equations that capture some aspect of a more compli- 
cated situation. In practice, some of the equations are exact (such 
as for example the mass and energy balance equation). On the 
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other hand, some of the physical variables are volume or mass 
averaged (e.g., the mean gas and temperature of the interstellar 
medium, respectively), or finally computed at some fiducial ra- 
dius of the assumed gas distribution. From this point of view, the 
specific galaxy and dark matter halo profiles do not enter directly 
into the code, but are needed to obtain realistic mean values to 
be used in the equations. 

A preliminary investigation of a physically based one-zone 
toy model has been done in SOCS, where it is shown that for a 
typical quasar spectral energy distribution (SED) the final Mbh 
produced by feedback clearly reproduces the observed Mbh-0" 
relation. The new models discussed here contain important im- 
provements with respect to SOCS. First of all, the modelization 
of accretion onto SMBH. In this paper we examine how the ef- 
fects of radiation pressure due to the Thomson electron scatter- 
ing modify in a self-consistent way the spherical Bondi flow. 
Another improvement is the treatment of Type la supernova rates 
and mass losses due to the evolution of the galactic stellar popu- 
lation, as we now adopt the Kroupa (2005) initial mas s function 
(IMF) coupled with the evolutionary prescription of [Marastonl 
(2005 1). Finally, th e dark matter halo is now described by a finite- 
mass Tjaffd (Il983h distribution instead of a singular isothermal 
sphere. 

We focus on a scenario in which the masses of the central 
SMBH and the host galaxy grow in a dark matter halo, which is 
replenished by accretion of gas of cosmological origin. We fol- 
low star formation and we also consider the mass return from 
the evolving stellar populations. The combined effect of SNIa 
heating and radiative feedback, during episodes when the lumi- 
nosity from the central SMBH approaches its Eddington limit, 
heats and drives much of the remaining gas out of the galaxy, 
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limiting both future growth of the SMBH and future star forma- 
tion to low levels. We do not consider the merging phenomenon, 
and we restrict ourselves to the evolution of an isolated galaxy: it 
is well known that significant AGN activity may also be present 
in isolated systems, because of the large amounts of gas pro 
duced by passively evolving stellar populations (e.g., Mathews 



1983; Ciotti et alJl99lHBregman & Parriott 2009; K avirai et al 
2010: ICisternas et al.l 120101) . More specifically, the one-zone 
models discussed in this paper are developed as an alternative 
a pproach to studying the " feedback modulated accretion flows" 
o flCiotti et al.l (12009, 2010). There are two major differences be- 
tween the present approach and hydro-simulations. The first, and 
most obvious, is the impossibility to resolve physical phenom- 
ena associated with specific length and timescales. The second is 
that we consider heating feedback onl y, while in the current ver - 
sion of the hydrodynamical code bv ICiotti et al.l (120091 [2010), 
both mechanical feedback, and radiation pressure are also con- 
sidered. However, despite these shortcomings, with the present 
models we can attempt to simulate the process of galaxy for- 
mation, which is beyond the current possibilities of the hydro- 
simulations and at the same time we can explore the parameter 
space. 

The evolutionary scenario that we consider here addresses 
several key observational findings. First, that giant ellipticals 
are old - they end their period of vigorous star formation 
early in cosmic time, since the radiative output from the cen- 
tral SMBHs limits (in cooperation with the energetic input due 
to star formation) the gas content to be at levels for which on- 
going star formation is minimal. Secondly, gas-rich, actively 
star-forming galaxies at redshift z ~ 3, including Lyman break 
galaxies and bright sub millimeter SCUBA galaxies, generally 
exhibit AGN act ivity (Steidel et alJl2002t lAlexander et al 1120031: 
lAlexanderll2009t iLehmer et al.ll2004t iDonlev et al.ll2010l) . indi- 



cating that their central SMBHs continue to grow. This sug- 
gests that the formation of a spheroid probably closely pre- 
ceeds a quasar shining phase, as verified by spectroscopic obser- 
vations indi cating that quasars occupy metal-enriched environ- 
ments (e.g ., Hama nn & Ferla nd 1999. ISchawinski et al.l [20091 
I Wild et all l2010r The redshift evolution of the quasar emis- 
sivity and the star formation history of spheroids is thus ex- 
pected to have evolved roughly in parallel since z ~ 3, which 
is also consistent w ith observations (e.g., Hai man et al.l 12004, 
iHeckman et al.ll2004T) . Among the most important observational 
predictions of the mod el is the length o f the so-called "obscured 
accretion phase" (e.g.. lComas tri 2004). This phase is defined as 
the period of time when a high column density is associated with 
a high accretion rate onto the central SMBH. We also study the 
relation between the duration of the obscured phase and the cor- 
responding "cold phase" (defined by a low mass-weighted gas 
temperature), and how they depend on the adopted parameters. 
Finally, in a large set of models we explore how the final SMBH 
mass is related to the final stellar mass, as a function of the dark 
matter halo mass, the amount of cosmological infall and the ac- 
cretion rate of the gas, the SMBH accretion radiative efficiency, 
and finally the initial SMBH mass. We find that the models are 
in close agreement with the observed Magorrian relation when 
we assume high efficiencies (e ~ 0.1) and relatively low initial 
SMBH masses (M BH < 1O 6 M ). 

This paper is organized as follows. In Sect. 2, we present 
the physics adopted in the simulations, with special emphasis on 
the differences from and improvements on SOCS. Section 3 is 
devoted to the description of the adopted self-regulated accretion 
model, while in Sect. 4 we present our findings. The main results 
are summarized and discussed in Sect. [5] while several useful 



interpolating functions that we use to compute the stellar mass 
losses and the SNe la rate are presented in the Appendix. 

2. The model 

The galaxy model and the input physics adopted for the simu- 
lations have been improved with respect to SOCS, in particular 
the treatment of SNIa heating, the mass return rate from evolving 
stars, and the accretion onto the central SMBH. Several aspects 
of this model were described in SOCS and Ciotti & Ostriker 
(2007), and we present here only the relevant modifications. 



2. 1 . The unchanged physics 

For completeness, we briefly summarize the aspects of the in- 
put physics that remain unchanged with respect to SOCS. The 
differential equation for the gas mass balance is 



M g = M M - AT. - / Edd M BH - M e: 



(1) 



The first source term on the right-hand side (r.h.s.) describes 
the cosmological infall in a pre-existing (and time-independent) 
dark matter halo 

Tinf 



Mm = 



(2) 



where M- m f is the total gas mass accreted during the simula- 
tion (but in general not equal to the final stellar galaxy mass 
M„), and tm is the characteristic infall timescale. A more ap- 
propriate description of the cosmological gas infall could be ob- 
tained by multiplying the r.h.s. of Eq. (f2]i by a factor t/Twi (e.g., 
iJohansson et al.ll2009l) . so that the total infalling mass is unaf- 
fected, while the infall starts in a more gentle way. For com- 
pleteness, we performed some simulation with this different de- 
scription, and we found (as expected) that for reasonable values 
of Tj n f the results are not much affected. Therefore, for consis- 
tency with SOCS we retain Eq. 0. The second source term is 
the net star-formation rate 



m, = mi - m:. 



where 



MI = 



(3) 
(4) 



is the instantaneous star-formation rate, and M™ is the mass re- 
turn by the evolving stellar population (Appendix IA.U ; follow- 
ing SOCS, in the simulations a„ is fixed to be 0.3. The charac- 
teristic star-formation timescale t» is defined as 



t, = max(T dyn , Toooi), 



(5) 



where the dynamical timescale Td yn is given by Eq. (fT8l ). and the 
mean gas cooling time t coo i is estimated to be 



^"cool 



Ec 



(6) 



In the equation above, Ec is the fiducial cooling rate given by 

* *m /Pgf X(l+X) t ,^ _ 

E c = n e n ? A(T) = — A(T), (7) 

\m p ) 2 

where p g is the instantaneous mean gas density (see Sect. 12. 31 ). 
A ^ = ^^IT 18 + 2 - 706 x l(T 47 r 2 - 976 erg s-'cm (8) 
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is the cooling functio n (Mathews & Breg manl 119781 : see also 
ICiotti & Ostrikerl20 1 ), and X = 0.7 is the hydrogen mass abun- 
dance (for simplicity we assume complete ionization). In agree- 
ment with Eqs. © and (|7), the mean gas internal energy is 



E = 



3p g k B T 
2fim p 



(9) 



where fi = (0.25 + 1.5X + 0.251T 1 - 0.62 is the mean atomic 
weight. It follows that 



^"cool 



fiX(X+l)M g A{T) 



(10) 



where M g is the instantaneous value of the total gas mass and r g 
is the gas distribution scale radius (Sect. 12. 3b . The third source 
term is the total accretion rate onto the SMBH 



(ID 



The first term describes gaseous accretion (Sect. [3), whereas the 
second term represents the contribution by the coalescence of 
stellar remnants of massive stars, as discussed in SOCS. The nu- 
merical coefficient /ndd ~ 10~ 2 - 10~ 3 needs to be implemented 
in the one-zone code (which by definition is unable to model 
the different spatial scales of the problem) to represent the ob- 
served time variation of quasars. In practice, /edd represents the 
"duty -cycle", and its value is constrained by both observations 
(e.g.lYu & Trem aine 2002, Haiman et al. 2004) and simulations 
(ICiotti & Ostrikerl2007LICiotti et al.l2009ll2010l) . When the ther- 
mal energy of the interstellar medium (ISM) of the galaxy is high 
enough, the gas is able to escape from the dark matter potential 
well, at a fiducial escape rate computed as 



M„ 







T ^ ^esc^v 



(12) 



The parameter 77 esc is of the order of unity, while the expression 
for r v ; r is given in the next section. Finally, the escape character- 
istic timescale is 

2r„ 

Tesc = — , (13) 

where c s is the speed of sound and r g is the scale length of the 
gas distribution (Sect. 12.31 ). Energy input into the ISM is pro- 
vided by the thermalization of supernova ejecta (both SNII and 
SNIa). The treatment of SNII is the same as in SOCS, and the 
new numerical treatment of SNIa is described in Appendix lA.2l 
Additional contributions to the ISM energetics come from the 
thermalization of red giant winds and radiative feedback due to 
accretion onto the SMBH. As in SOCS, we adopt the average 
quasar SED obtained from the CRB supplemented by informa- 
tion from individual objects. We recall that the UV and high en- 
ergy radiation from a typical quasar can photoionize and heat 
a low density gas up to an equilibrium Compton temperature 
(7c ~ 2 x 10 7 K) that exceeds the virial temperatures of giant 
ellipticals. Following SOCS, we also consider adiabatic cool- 
ing in the case of gas escaping and heating/cooling due to in- 
flow/outflowing galactic gas. The gas temperature is therefore 
determined at each time-step by integrating the equation of the 
internal energy per unit volume 



E - -Eh.SN + £h,w + £h,AGN - Ec+ £ad + 3 



Minf^L - M ( 



esc^s 



(14) 



where £h,sn is the energy due to SNIa and SNII, £h,w describes 
the thermalization of red giant winds, £h,agn is the AGN heat- 
ing, £ a d is the adiabatic cooling in the case of galactic winds, 
and A is a dimensionless parameter (0.25 < X < 1). Finally, we 
force the gas to remain above 10 4 K. 

2.2. The dark matter halo 

For simplicity, the dark matter (DM) halo in our model is the 
only contributor to the gravity of the galaxy. In addition, the gas 
and the stellar density distributions in the simulations are as- 
sumed to be proportional to the local (unevolving) DM density 
distribution, modeled as a Jaffe (1983) profile (Sect. l231 > 



Ph(r) = 



M h 



4n r 2 (r h + r) 2 ' 



(15) 



For a DM halo total mass M h, the scale length rj, is fixed fol- 
lowing iLanzonTetal] d2004l) . In that paper, 13 massive DM ha- 
los obtained from high-resolution cosmological N-body simu- 
lations were carefully analyzed. In particular, for each halo the 
"overdensity radius" r& (sometimes called virial radius) was de- 
termined, and it was shown that r& does not differ significantly 
from the true virial radius r v ; r of the system (with a scatter < 
20%). Averaging oy er the 13 clusters reported in Table 2 in 
lLanzoni et al.l d2004l) . we found that 



M h a 0.03>|, 



(16) 



where Mh is in 1O 9 M and is in kpc units. For the Jaffe profile 
r v ; r = 2rh, and from Eq. ( [ToT i we can link the scale-length to 
Mh by assuming r wu = r&. We finally recall that the virial (3D) 
velocity dispersion associated with Eq. ( fT3l ), is given by 



2 GM h 



2r h 



(17) 



The halo mean circular velocity is estimated as v\ - 2cr 2 ^, and 
the dynamical time required in Eq. (0 is therefore given by 



Tdyn — 



2nr h 



(18) 



2.3. The galaxy model 



One of the new aspects of the model evolution investigated in 
this paper is the time extent of the Compton thin and Compton- 
thick phases, as defined at the end of this section. In SOCS, the 
adopted gas density distribution was a singular isothermal sphere 
and the major drawback of this choice for the evaluation of aver- 
age quantities is the divergence of the total mass. For this reason, 
we now use the more realistic Jaffe density profile, i.e., we as- 
sume that at each time the gas density distribution is 



Pg( r ) = 



An r 2 (r a + r) 2 ' 



(19) 



where r g is the gas scale radius and M g is the instantaneous 
value of the total gas mass, obtained by integrating Eq. (Q}. For 
simplicity, r g - during the whole simulation. However, it is 
possible to generalize the present approach and consider a time- 
dependent r g . 

The mean gas density, which is needed in Eqs. (0 and (0, 
is evaluated to be the instantaneous mean value within the half- 
mass radius r„ 

3M„ 

P g = — f , (20) 
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while the virial gas temperature r v i, is obtained from the hydro- 
static equilibrium and the Jeans equations solved in the gravita- 
tional potential of the DM halo neglecting the self-gravity of the 
gas. Simple algebra then shows that 



fim p GMh 
6k B r„ 



(21) 



To quantify the relative importance of the optically thin and 
thick phases, the code computes the gas column density as 



M Pg (R) 
2nR 2 ' 



(22) 



where Mp g (R) is the projected gas mass enclosed in a circle of 
radius R, and the factor of 2 at the denominator takes into ac- 
count that only one side of the gas column actually obscures the 
center. Therefore, the fiducial column density depends not only 
on the total gas mass, but also on t he aperture radius adopted. 
To simulate observational work (e.g. lRisaliti et al.lll999h . we de- 
cided to define R to be one hundredth of the effective radius R^\, 
i.e., in Eq. fl22]> R = 0.0074r g . From Eq. dATTb 



(N H ) = 32.7- 



(23) 




and r g is fixed during the simulation, (N H ) depends only on M g . J^g- l^he self " re g ulated Bondi accretion rate M* obtained from 



We refer to the "thick" phase if (Nu) > 10 cm" , while galax- 
ies with column density (Nh) > 10 22 cm" 2 are considered ob- 
scured. 



Eq. (IBTt is represented by the heavy solid line. The dotted line is 
the accretion rate determined by Eq. d251 l. Note how for very low 
and very high classical accretion rates M B , the two prescriptions 
coincide. 



2.4. The code 

The simulations are performed by numerical integration of the 
previous evolutionary differential equations using a forward 
Euler scheme. The numerical code is based on the code devel- 
oped in SOCS. The time step is defined as the minimum among 
the characteristic times associated with the different physical 
processes 

(24) 



St = B min (f inf , f„, f g , t BH , t E ) , 



where the different subscripts indicate the specific aspect of the 
physics involved and, in general, for a quantity X[, ti = Xj/\Xi\. 
The dimensionless coefficient B < 1 is used to improve the code 
accuracy. We performed several test simulations and we found 
that B = 0.1 leads to rapid convergence and excellent agree- 
ment with results obtained in SOCS, when used with their model 
galaxies and input physics. 

3. Self-regulated Bondi accretion 

As often assumed in the literature, the SMBH accretion rate is 
determined as the minimum between the Bondi accretion rate 
Mg and Eddington limit M Edd 



M BH ,acc = min(M Edd , M B ). 
In the previous formula, the Eddington accretion rate is 

^Edd 



M Edd = 



where 0.001 < e < 0.1 is the accretion efficiency, and 

AncGM-snnm-p 



^Edd 



(25) 



(26) 



(27) 



For the Jaffe model, R c as 0.74r„ 



is the Eddington luminosity and cr T is the Thomson cross- 
section. The classical Bondi accretion rate is 



A c 4tt G 2 M 2 h pooC s i 



(28) 



where and c s oo are the gas density and the speed of sound 
at infinity, respectively, and A c is a d i mensionless c oefficient of 
the order of unity (e.g., lBondilll952l lKroliklll999l) . In eq (f28]i, 
radiative effects are not taken into account, so that in Eq. ( FSB! 
there is a sharp transition between the pure hydrodynamical 
and radiation-dominated regimes. Fortunately, in the optical thin 
regime dominated by electron scattering, it is possible to extend 
the classical Bondi accretion solution and take into account the 
radiation pressure effects, so that the transition between Bondi- 
limited and Eddington-limited accretion can be described in a 
more consistent way. After solving this pro blem, we disco vered 
that is had alr eady been fully described bv lTaam et al.l(l 19911) and 
iFukud (120011) . and here we just summarize the main features of 
the modified accretion model. The result is obtained basically, 
by noting that in the optically thin regime the radiation pressure 
scales as 1 / r 2 , thus reducing the effective gravitational force at 
each radius by the same amount. By imposing self-consistentcy, 
i.e. by requiring that the effective accretion rate Mg is deter- 
mined by the effective gravity, one finally obtains the following 
expression for the modified coefficient 



*c,eff 



1 - 



M e 

M E dd) 



(29) 



Therefore, the self-consistent Bondi accretion rate satisfies the 
equation 



M e R =M B 



1 



M Eddj 



(30) 
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Fig. 2. The time evolution of relevant quantities of the RM. Panel 
a: gas temperature; the virial temperature of the galaxy is the 
horizontal dashed line. Panel b: bolometric accretion luminos- 
ity. Panel c: gas density at the Bondi radius (black line) and the 
mean gas density (green dotted line). Panel d: mean gas column 
density (Nh) at aperture radius of 0.088 kpc. 

that can be solved for M B . After discarding the unphysical solu- 
tion, we have 

• - 1[. - n ~ 2 i _ M Edd 

m = —. = - 2 + r - V 4r + r 2 , r = — : — . (31) 

M Edd 2L J M B 

For r — > (high accretion rates) m — > 1, so that M B tends to 
the Eddington accretion, while for r — > oo (low accretion rates) 
m — > 0. The full solution is indicated in Fig.[T]by the heavy solid 
line: we note how Eq. d25l > overestimates accretion onto SMBH 
in the range 0.1 < M B /M Edd < 100. 

4. Results 

In the following sections, we present the main properties of a set 
of simulations, focusing on three important issues and compar- 
ing the new results with those in SOCS. However, before dis- 
cussing the whole set of the new models, in the next section we 
present in detail the evolution of a representative model and three 
of its possible variants. 

4.1. The evolution of the reference model and some of its 
variants 

Figure|2]shows the time evolution of important quantities of our 
galaxy reference model (RM). This model is characterized by a 
DM halo of total mass Mh = 4 x 10 11 M Q , which corresponds to 
a halo scale length rj, = 11.86 kpc and R e = 8.83 kpc, a fidu- 
cial circular velocity calculated according to Eq. dTTT > of ~ 270 
km/sec, and a characteristic infall time 2 Gyr. The assumed total 
mass of the cosmological gas infall is M m f = 10 11 M , corre- 
sponding to a dark-to-total mass ratio of 80% (provided that all 




5 10 15 

t/Gyr 



Fig. 3. The SMBH accretion history of the RM. Top panel: time 
evolution of the SMBH mass as computed according to SOCS 
(green line) and with the Bondi-modified accretion (black line). 
Bottom panel: SMBH accretion rate from SOCS (green line) and 
from our work (black line). The red solid line is the Eddington 
accretion rate for the SOCS model. 



the infalling gas forms stars). Other simulation parameters are 
a. = 0.3,/3bh,. = 1.5 x 10~ 4 , e = 0.1, tj sn = 0.8 5, = 2, and 
a stellar mass-to-light ratio of 5 (see Eq. IA.11I ). As in SOCS, 
the initial SMBH mass is 10 8 M Q ; the duty circle, / Edd , is fixed 
at 0.005. 

As can be seen from the comparison of Fig. [2] with Figs. 7 
and 8 in SOCS, the global evolution of the new models is qual- 
itatively very similar to those in SOCS. In particular, the time 
evolution of the model, from the beginning up to ~ 9 Gyrs, is 
characterized by a cold phase (we assume that a cold phase cor- 
responds to Tgas < 10 5 K) of high density and low temperature. 
The gas density at the Bondi radius remains between 1 and 10 2 
particles per cm 3 , while the mean gas density is ~ 10 _1 - 10~ 2 
particles per cm 3 . At the beginning of the cold phase, about 2.17 
Gyrs are spent in the Compton-thick phase (see Fig . [2]d) . The re- 
maining part of the cold phase is obscured, while the accretion 
luminosity remains high, at Lbh <; 10 46 erg/sec. As soon as the 
mean gas density decreases (Fig. [2}; green dotted line) the cool- 
ing becomes inefficient, the gas heating dominates, and the tem- 
perature increases to the virial temperature. The total duration 
of the cold phase is ~ 9 Gyr, which corresponds to a decrease 
in (Nh) of about 2 orders of magnitude. The obscured phase 
lasts ~ 9.6 Gyr and the average column density is then around 
(Nh) ~ 10 20 crrT 2 . The durations of the cold and the obscured 
phases are very similar, and indeed, the two phases are related. 
Until the gas density is high, the gas can radiate efficiently the 
energy input due to the AGN, and its temperature remains below 
10 5 K. Because of this large amount of cold gas, the star forma- 
tion proceeds at high rates, until most of the gas is consumed. 
At this point, the radiative cooling of the gas becomes ineffi- 
cient and the heating due to the AGN causes the gas temperature 
to increase. The combined effects of the gas consumption and 
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Table 1. Final properties of the RM and its variants at 15 Gyrs. 



Model 


M, 


M csc 




M BH /M, 




10g(iVH> 




cold 


AfcT 


Af b s 




RM 


93.22 


6.56 


0.256 


2.75 x 10" 


3 


20.05 


10.42 


7.94 


2.17 


9.57 


3.26 


RM, 


92.47 


6.62 


0.945 


1.02 x 10" 


-2 


20.04 


10.37 


7.91 


2.13 


9.56 


3.31 


RM 2 


42.65 


7.21 


0.214 


5.01 x 10" 


-3 


19.63 


8.42 


5.78 


0.00 


10.33 


4.67 


RM 3 


192.97 


6.69 


0.321 


1.67 x 10" 


-3 


20.46 


12.04 


9.69 


4.92 


8.17 


1.91 


RM 4 


90.59 


5.55 


0.288 


3.18 x 10" 


-3 


22.81 


14.20 


10.47 


0.00 


14.99 


0.01 



Note - All masses are in 10 9 M Q units. Mbh is the final black hole mass; M, is the final stellar mass; M esc is the total escaped gas mass. The average 
column density is in cm" 2 . Af CT > Af obs , and Af unobs represent the durations (in Gyr) of the Compton-thick phase, of the obscured phase, and the 
unobscured phase, respectively. Finally, AfJ.. and Af 2 o]d are the durations (in Gyr) measured assuming threshold temperatures of 10 5 K and 5 x 10 4 
K, respectively. 



the increase temperature stop the star formation. A specific fea- 
ture of the cold phase can be seen in Fig. [2] where characteristic 
temperature oscillations are apparent. These oscillations, already 
presented and discussed in SOCS, are due to the combined effect 
of gas cooling and AGN feedback. They terminate when the gas 
density falls below some threshold determined by the cooling 
function and because less and less gas is produced by stars and 
accreted by the galaxy, while SNIa heating declines less strongly. 
When the density is high (at early times), the cooling time is in- 
stead, very short, and AGN heating is radiated efficiently. The 
two competitive effects produce the temperature (and density) 
oscillations. In the hydro-simulations, the spatial and temporal 
structures of these oscillations is quite complicated, as feedback 
and cooling act on several different spatial and temporal scales 
(from a month to 10 Myr). In the one-zone models, the oscilla- 
tions are instead dependent on the "duty-cycle" parameter fm& 
in Eq. (Q]), which is constrained by observational and theoretical 
studies (see Sect. 12. U . 

To summarize, the observational properties of the RM 
would correspond to a system initially Compton-thick, that then 
switches to an obscured phase, and in the past 3.2 Gyrs the 
galaxy is unobscured. The SMBH accretion history is shown in 
Fig. [3] where it is compared to that predicted by the SOCS for- 
mulation in a model otherwise identical to the RM. The major 
difference between the two models is in the final value of Mbh: 
in particular, Eq. ( |25*1 ) would predict a present-day Mbh a factor 
~ 1 .5 higher than that for modified accretion. This is because the 
accretion in SOCS remains for almost all of the galaxy life, at the 
Eddington limit, while in the new treatment the self-regulation 
maintains the accretion at lower rates. 

Before presenting the results of a global exploration of the 
parameter space, we focus on a few obvious questions. For ex- 
ample, what happens if we keep all the model parameters fixed 
and reduce the radiative accretion efficiency? What happens if 
we increase or decrease the total mass of infalling gas? Or if we 
double the gas infalling time? Of course, these simple examples 
do not cover all the possible cases. However, these models will 
provide a guide as we investigate the parameter space. The rele- 
vant properties of the additional "reference" models are listed in 
Tabled 

The first variant of the RM model is model RMi obtained 
by reducing the radiative accretion efficiency from e = 0.1 to 
e = 0.001. Overall, this reduced efficiency model passes through 
the same evolutionary phases as model RM: an initial Compton- 
thick phase, followed by an obscured phase, and finally a low- 
density unobscured phase. The initial cold high density phase is 
also very similar to that of RM in Fig. [2] The most important 
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Fig. 4. Evolution of the SMBH mass and accretion rate com- 
puted according to Eq. d25l > (green line, SOCS) and Eq. d30b 
(black line) for a model with M inf = 1.25 x 10 I0 M o , M h = 
5 x 1O 1O M , and initial SMBH mass of 10 8 M o . From left to right; 
e = 0.1, e = 0.01, and e = 0.001. The red line represents the 
Eddington limit for the SOCS model. 

and expected difference is in the final mass of SMBH, which 
is higher by a factor of a: 3.5 than the RM. A reduction in the 
radiative efficiency also produces a slightly larger amount of es- 
caped gas, because of the shorter cold phase. As a consequence, 
a larger amount of gas produced by the evolving stars is lost from 
the galaxy during the low-density, hot phase. 

The effect of a reduction in M;„f is explored in model 
RM2, while all the others parameters are the same as in RM. 
Qualitatively, the evolution is again very similar to that of RM in 
Fig. [2] The main difference is the absence of the initial Compton- 
thick phase because of the lower gas density. Both the final stel- 
lar mass and final SMBH mass are lower than in model RM, as 
can be seen from Table Q] Overall, model RM2 does not display 
remarkable or unexpected properties. The only noticeable as- 
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pect is that the escaped mass is higher than in model RM, but 
this is again due to the shorter cold phase. 

A complementary model to RM2 is RM3, where we dou- 
ble the value of M; n f while maintaining all the other parame- 
ters identical to those of model RM. The qualitative evolution is 
again similar to that of RM in Fig.[2l but, at variance with RM2, 
the Compton-thick phase is now present. Not surprisingly, the 
Compton-thick phase in model RM3 extends for a longer period 
(~ 5 Gyrs) than in the other variants of RM (see Table Q]), and 
its total mass of new stars is also the highest. However, the total 
mass that escapes is not as high as one would expect, the infall 
mass being a factor of 2 higher than in the RM yet the escaped 
mass being almost the same. The main reason for this behaviour 
is the very massive star formation, which is almost double that of 
RM. We note that model RM 3 is the closest, in the RM family, 
to the galaxy models studie d in the hydrodynamical simulations 
of ICiotti et al.1 (120091 1201 Oh . as far as the final stellar mass and 
the final Mbh are concerned. 

We finally discuss model RM4, which is identical to model 
RM, but has twice as long an infall time. The main effects are the 
longest cold phase in RM family (see Table [TJ, and the absence 
of the initial Compton-thick phase. This latter characteristic is 
due to the time dilution of the infalling gas density, which pre- 
vents the possibility of reaching very high column density val- 
ues. Overall, the final SMBH mass is not affected significantly 
by the extended infall phase, but its escaped gas masses is quite 
high because the total mass in new stars is (as expected) lower 
than in the RM model. 

To summarize, these preliminary experiments have revealed 
that, at fixed dark matter halo, sensible variations in the input 
parameters do not produce remarkable different results. The rel- 
evant differences are found mainly in the amount of star forma- 
tion and the different durations of the cold and obscured phases, 
while the final SMBH mass appears to be mainly affected by the 
value of the radiative efficiency e. 

4.2. Exploring the parameter space 

We now present the general results (summarized in Tables|2]and 
of our examination of the parameter space. The SMBH ini- 
tial mass in the models in Table [2] is 1O 8 M , while in Table [3] it 
is 1O 5 M . From the astrophysical point of view, the first choice 
mimics a scenario in which the central SMBHs are already quite 
massive at the epoch of galaxy formation, while in the second 
case the main growth occurs with galaxy formation. Each of the 
5 main families of models (Mi, M5) consists of galaxy mod- 
els characterized by the same dark halo mass, Mh, ranging from 
10 10 M o (the Mi family) up to 10 12 M o (the M 5 family). In each 
family, we have investigated 6 models with different values of 
infalling gas-to-dark matter ratio a = Mj n f /Mh, and finally dif- 
ferent radiative efficiencies e in the self -regulated Bondi accre- 
tion. In practice, in each family we explore the effects of differ- 
ent total infalling gas mass and different efficiencies (spanning 
the commonly accepted range of values) at fixed Mh. 

Almost independently of Mh and the initial value of Mbh, we 
can recognize some general trends. For example in each group 
of 3 models characterized by identical parameters but decreasing 
e, it is apparent how the total amount of stars formed decreases 
at decreasing e, while the final Mbh increases. In the total mass 
budget of the galaxy, an important quantity is the total mass of 
gas ejected, M esc . As expected, we find that massive galaxies, 
with larger infall gas masses also eject more mass. However, in 
all cases, the escaped gas mass is much lower than the final stel- 
lar mass, i.e., ~ 10% of it. 



Finally, higher final gas masses are obtained, at fixed Mh, for 
higher M; n f and lower efficiencies, because of the less effective 
feedback. We note however that the product of e and the mass 
accreted by the SMBH decreases with decreasing e, i.e., even 
though the accreted mass is higher, the integrated energy output 
is lower, and so is the feedback effecfl 

A qualitative illustration of the effect of the reduction in e 
on the SMBH accretion history is shown in Fig. [4] When the ef- 
ficiency decreases, the final mass of the SMBH increases and, 
we indicate by the green line, the corresponding evolution of 
identical models but with the SOCS treatment, i.e., where the 
accretion is determined by the minimum of Mb and MecM- As ex- 
pected, the differences in Mbh reduce with decreasing e, because 
the Eddington accretion regime in the SOCS models (when the 
major differences from the Bondi-modified case are established), 
becomes less and less important, and M B approches Mb- 

All models have a transition phase from obscured to unob- 
scured. The Compton-thick phase is present in almost all models 
that have a cold phase (18 of 30 simulations); however, in the 
less massive set of models (Mi) the Compton-thick phase is ab- 
sent. This is consistent with the scenario in which the majority 
of the most massive black holes spend a sig nificant amount of 
time growing in an earlier obscured phase (see lKellv et all 2010: 
iTreister et al]|2010l) . 

Finally, the mean SMBH-to-sta r ratio is not far fro m the 
value inferred from the present-day Mago rrian et al.l dl998l) re- 
lation, though on the high mass side. This is not surprising, be- 
cause of our use of a one-zone model. In any case, we empha- 
size that the SMBH feedback in the models was able to remove 
(in combination with star formation and/or gas escape) most of 
the infalling gas (which, if accreted onto the central SMBH in a 
cooling-flow like solution, would lead to the final SMBH masses 
being ~ 2 orders of magnitude higher than the observed ones). 

Independent of the particular characteristics of the single 
models, Fig. [5] clearly illustrates that the final SMBH masses, 
in models with a quite high initial Mbh (Table 0, are higher 
than implied by the observed Magorrian relation. For this rea- 
son, we explored other families of models, in which the initial 
mass of the SMBH was reduced to 1O 6 M , 10 5 M o , and 1O 3 M . 
From Fig. [5] it is apparent how galaxy models with initial SMBH 
masses < 1O 6 M (circles) closely agree to within l<x of the ob- 
served dispersion with observations, but only when the radiative 
efficiency is high, i.e. e = 0.1. It is even more remarkable that 
the final SMBH mass is proportional to the final M„, as all the 
models started with the same initial SMBH mass. This strongly 
indicates that co-evolution is the most plausible explanation of 
the proportionality between M„ and Mbh, ip line with other ob- 
servational evidence (e.g., see lHaiman et al . 200-1 and references 
therein). 

To compare our present findings with those of a complemen- 
tary hydrodynamical approach, in Fig. [5] we plot the results of 
the high-resolution hydro-simulations developed by Ciotti et al. 
(2010, Table 1, Cols. 5 and 6) for a representative galaxy with 
stellar mass of ~ 3 x 1O U M . The final Mbh are represented 
by crosses, and the different (luminosity weighted) radiative ef- 
ficiencies are in the range 0.003 < e < 0.133. We note that 
the crosses indicate almost all the hydrodynamical models that 
have been studied in detail so far, and this shows the importance 
of one-zone models as a complementary approach to exploring 
the parameter space. We also note that the final mass of the 



2 The higher SMBH masses account for the slightly lower final stellar 
masses. The mass conservation of the code is quite remarkable, consid- 
ering the amount of input physics involved (M m f = M, + M esc + Mbh)- 
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Table 2. Final properties of models with self-regulated Bondi accretion with an initial SMBH mass of 10 8 M o . All masses are in 
units of 1O 9 M . 





a 


6 


M, 


M esc 


M B h 








log<iV H > 


Af cold 




Af CT 


Afobs 
















IV 


1 n 
l = 10 
















2.5 


0.25 


0.1 


2.40 


0.07 


0.132 




5.51 x 10" 


-2 


18.95 


10.74 


8.39 


0.00 


11.29 


3.71 


2.5 


0.25 


0.01 


2.38 


0.07 


0.154 




6.45 x 10" 


-2 


18.94 


10.69 


8.35 


0.00 


1 1.23 


3.77 


2.5 


0.25 


0.001 


2.37 


0.07 


0.158 




6.67 x 10" 


-2 


18.94 


10.68 


8.35 


0.00 


11.22 


3.78 


5.0 


0.5 


0.1 


4.88 


0.07 


0.144 




2.95 x 10" 


-2 


19.36 


12.18 


9.96 


0.00 


12.70 


2.30 


5.0 


0.5 


0.01 


4.85 


0.07 


0.177 




3.65 x 10" 


-2 


19.35 


12.13 


9.92 


0.00 


12.64 


2.36 


5.0 


0.5 


0.001 


4.84 


0.07 


0.184 




3.80 x 10" 


-2 


19.35 


12.12 


9.92 


0.00 


12.64 


2.36 














M 


— 

2 — jy) 
















12.5 


0.25 


0.1 


11.23 


1.20 


0.167 




1.49 x 10 


-2 


19.43 


10.51 


8.11 


0.00 


1 1.45 


3.55 


12.5 


0.25 


0.01 


11.12 


1.21 


0.257 




2.32 x 10" 


-2 


19.42 


10.45 


8.07 


0.00 


11.39 


3.61 


12.5 


0.25 


0.001 


1 1.07 


1.23 


0.289 




2.61 x 10" 


-2 


19.42 


10.45 


8.06 


0.00 


11.38 


3.62 


25.0 


0.5 


0.1 


24.10 


0.79 


0.190 




7.90 x 10 


-3 


19.84 


12.15 


9.87 


2.31 


10.55 


2.14 


25.0 


0.5 


0.01 


23.96 


0.80 


0.328 




1.37 x 10" 


-2 


19.83 


12.10 


9.84 


2.29 


10.51 


2.20 


25.0 


0.5 


0.001 


23.92 


0.80 


0.370 




1.55 x 10" 


-2 


19.83 


12.10 


9.84 


2.29 


10.52 


2.20 












M 


— 1UU 
















25.0 


0.25 


0.1 


23.07 


1.82 


0.189 




8. 18 x 10 


-3 


19.64 


10.53 


8.11 


0.00 


11.51 


3.49 


25.0 


0.25 


0.01 


22.91 


1.83 


0.346 




1.51 x 10" 


-2 


19.63 


10.48 


8.07 


0.00 


1 1.45 


3.55 


25.0 


0.25 


0.001 


22.84 


1.84 


0.414 




1.81 x 10" 


-2 


19.63 


10.48 


8.07 


0.00 


1 1.45 


3.55 


50.0 


0.5 


0.1 


48.34 


1.51 


0.220 




4.55 x 10" 


-3 


20.05 


12.14 


9.84 


3.33 


9.58 


2.08 


50.0 


0.5 


0.01 


48.10 


1.52 


0.452 




9.40 x 10 


-3 


20.04 


12.09 


9.80 


3.32 


9.55 


2.14 


50.0 


0.5 


0.001 


48.00 


1.53 


0.536 




1.12 x 10" 


-2 


20.04 


12.09 


9.81 


3.31 


9.55 


2.14 












M 


t — JUU 
















125.0 


0.25 


0.1 


115.84 


8.91 


0.273 




2.36 x 10 


-3 


20.1 1 


10.37 


7.87 


2.51 


9.33 


3.15 


125.0 


0.25 


0.01 


115.38 


8.87 


0.776 




6.73 x 10" 


-3 


20.10 


10.32 


7.84 


2.50 


9.29 


3.21 


125.0 


0.25 


0.001 


114.95 


8.97 


1.104 




9. 61 x 10" 


-3 


20.10 


10.32 


7.84 


2.48 


9.31 


3.21 


250.0 


0.5 


0.1 


240.57 


9.03 


0.349 




1.45 x 10" 


-3 


20.52 


12.01 


9.65 


5.14 


8.01 


1.84 


250.0 


0.5 


0.01 


239.84 


9.06 


1.048 




4.37 x 10" 


-3 


20.51 


11.98 


9.63 


5.13 


7.98 


1.89 


250.0 


0.5 


0.001 


239.42 


9.12 


1.407 




5.88 x 10" 


-3 


20.52 


11.98 


9.63 


5.12 


7.99 


1.89 












M 5 


= 1000 
















250.0 


0.25 


0.1 


222.85 


26.74 


0.350 




1.57 x 10" 


-3 


20.31 


10.09 


7.53 


3.31 


9.18 


2.50 


250.0 


0.25 


0.01 


222.06 


26.68 


1.216 




5.48 x 10" 


-3 


20.30 


10.05 


7.51 


3.30 


9.08 


2.61 


250.0 


0.25 


0.001 


220.80 


27.16 


1.994 




9.03 x 10" 


-3 


20.30 


10.04 


7.50 


3.28 


9.10 


2.62 


500.0 


0.5 


0.1 


471.83 


27.47 


0.476 




1.01 x 10" 


-3 


20.72 


11.85 


9.44 


5.78 


7.77 


1.45 


500.0 


0.5 


0.01 


470.52 


27.66 


1.613 




3.43 x 10" 


-3 


20.71 


11.82 


9.43 


5.77 


7.71 


1.52 


500.0 


0.5 


0.001 


469.62 


27.85 


2.336 




4.97 x 10" 


-3 


20.71 


11.83 


9.43 


5.77 


7.72 


1.52 



Note - The halo mass M h is (10'°,5 x 10 10 , 10", 5 x 10", 10 12 ) M Q for models M[,M 2 ,M3,M 4 , and M 5 , respectively. The total infall gas mass 
is M- mf = aM h . M Bli is the final SMBH mass (also with the contribution of the initial SMBH mass). M, is the final stellar mass. M csc is the total 
escaped gas mass. The average column density is in cm" 2 . A?ct, A? b S , and A? uno b s represent the duration (in Gyrs) of the Compton-thick phase, 
the obscured phase, and the unobscured phase, respectively. Finally, A/' old and Af 2 o]d are the durations (in Gyr) measured assuming threshold 
temperatures of 10 5 K and 5 x 10 4 K, respectively. 



SMBHs, as computed in the hydrodynamical simulations, lead 
to an accurate reprodution of the Magorian relation but the com- 
parison with the one-zone models is delicate. In the hydrody- 
namical simulations (at variance with the one-zone models), the 
initial phases of galaxy formation are not simulated, and the fo- 
cus is on the maintaining low SMBH masses in the presence of 
stellar mass losses that, if accreted onto the central SMBH as in 
an undisturbed cooling flow, would produce a final BH mass of 
about a factor of ~100 higher than the observed ones. For this 
reason, in the hydrodynamical models the galaxy is already as- 
sumed to have formed, and the initial mass of the SMBH is just 
slightly lower than the mass predicted by the Magorrian relation. 



5. Conclusions 

This paper is a natural extension of a previous paper by Sazonov 
et al. (2005). From a technical point of view, several as- 
pects of the input physics have been improved. In particular, 
we have adopted a different description of the accretion rate, 
which now follows the self-con sistent modified Bondi theory 
(Taam et all 11991b iFukud 1200 ll) . instead of the minimum be- 
tween the Eddington and the (classical) Bondi rate. Moreover, 
the SNIa rate is now computed using a high-precision multi- 
exponential approximation for the time-kernel in the convolu- 
tion integral of the star formation rate, instead of the standard 
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Table 3. Final properties of models with self-regulated Bondi accretion with an initial SMBH mass of 10 5 M o . All masses are in 
units of 10 9 M o . 



M iaf 


a 


6 


M, 


M csc 




M BH /M, 




log<jV H > 


Af cold 




A/ct 


Af obs 




Mi = 10 


2.5 


0.25 


0.1 


2.43 


0.06 


0.002 


6.72 x 10- 


-4 


18.98 


10.94 


8.50 


0.00 


11.51 


3.49 


2.5 


0.25 


0.01 


2.41 


0.07 


0.026 


1.06 x 10- 


-2 


18.94 


10.70 


8.37 


0.00 


11.24 


3.76 


2.5 


0.25 


0.001 


2.36 


0.07 


0.071 


3.02 x 10- 


-2 


18.93 


10.64 


8.31 


0.00 


11.18 


3.82 


5.0 


0.5 


0.1 


4.93 


0.07 


0.003 


6.39 x 10- 


-4 


19.39 


12.36 


10.06 


0.00 


12.91 


2.09 


5.0 


0.5 


0.01 


4.88 


0.07 


0.044 


9.01 x 10- 


-3 


19.35 


12.11 


9.90 


0.00 


12.62 


2.38 


5.0 


0.5 


0.001 


4.81 


0.07 


0.113 


2.35 x 10- 


-2 


19.34 


12.09 


9.88 


0.00 


12.60 


2.40 


M 2 = 50 


12.5 


0.25 


0.1 


11.32 


1.17 


0.006 


5.51 x 10- 


-4 


19.45 


10.61 


8.16 


0.00 


11.55 


3.45 


12.5 


0.25 


0.01 


11.23 


1.19 


0.071 


6.36 x 10- 


-3 


19.43 


10.47 


8.09 


0.00 


11.39 


3.61 


12.5 


0.25 


0.001 


11.10 


1.20 


0.192 


1.73 x 10- 


-2 


19.42 


10.44 


8.06 


0.00 


11.37 


3.63 


25.0 


0.5 


0.1 


24.19 


0.78 


0.014 


5.60 x 10- 


-4 


19.86 


12.23 


9.92 


2.32 


10.63 


2.05 


25.0 


0.5 


0.01 


24.06 


0.79 


0.130 


5.41 x 10- 


-3 


19.83 


12.09 


9.84 


2.32 


10.48 


2.21 


25.0 


0.5 


0.001 


23.90 


0.79 


0.297 


1.24 x 10- 


-2 


19.83 


12.09 


9.83 


2.28 


10.51 


2.21 


M 3 = 100 


25.0 


0.25 


0.1 


23.18 


1.79 


0.012 


5.31 x 10- 


-4 


19.65 


10.60 


8.14 


0.00 


11.58 


3.42 


25.0 


0.25 


0.01 


23.06 


1.81 


0.119 


5.16 x 10- 


-3 


19.63 


10.49 


8.09 


0.00 


11.45 


3.55 


25.0 


0.25 


0.001 


22.87 


1.82 


0.300 


1.31 x 10" 


-2 


19.63 


10.47 


8.07 


0.00 


11.44 


3.56 


50.0 


0.5 


0.1 


48.44 


1.50 


0.026 


5.42 x 10" 


-4 


20.06 


12.19 


9.86 


3.34 


9.64 


2.02 


50.0 


0.5 


0.01 


48.24 


1.52 


0.212 


4.40 x 10" 


-3 


20.04 


12.08 


9.80 


3.33 


9.52 


2.15 


50.0 


0.5 


0.001 


48.00 


1.52 


0.452 


9.42 x 10" 


-3 


20.04 


12.08 


9.80 


3.31 


9.55 


2.14 


M 4 = 500 


125.0 


0.25 


0.1 


116.05 


8.82 


0.059 


5.07 x 10" 


-4 


20.12 


10.39 


7.88 


2.52 


9.35 


3.13 


125.0 


0.25 


0.01 


115.67 


8.85 


0.407 


3.52 x 10" 


-3 


20.10 


10.33 


7.85 


2.51 


9.28 


3.21 


125.0 


0.25 


0.001 


115.21 


8.86 


0.868 


7.54 x 10" 


-3 


20.10 


10.32 


7.84 


2.49 


9.29 


3.22 


250.0 


0.5 


0.1 


240.70 


9.02 


0.123 


5.10 x 10" 


-4 


20.53 


12.03 


9.66 


5.14 


8.03 


1.82 


250.0 


0.5 


0.01 


240.15 


9.03 


0.681 


2.84 x 10- 


-3 


20.51 


11.97 


9.63 


5.14 


7.97 


1.89 


250.0 


0.5 


0.001 


239.57 


9.06 


1.238 


5.17 x 10- 


-3 


20.51 


11.98 


9.63 


5.13 


7.98 


1.89 


M 5 = 1000 


250.0 


0.25 


0.1 


222.92 


26.81 


0.114 


5.12 x 10- 


-4 


20.31 


10.10 


7.54 


3.31 


9.21 


2.48 


250.0 


0.25 


0.01 


222.41 


26.71 


0.741 


3.33 x 10" 


-3 


20.30 


10.05 


7.51 


3.31 


9.08 


2.61 


250.0 


0.25 


0.001 


221.73 


26.64 


1.504 


6.78 x 10" 


-3 


20.30 


10.05 


7.51 


3.30 


9.07 


2.63 


500.0 


0.5 


0.1 


471.96 


27.48 


0.238 


5.05 x 10" 


-4 


20.72 


11.86 


9.45 


5.78 


7.78 


1.44 


500.0 


0.5 


0.01 


471.01 


27.55 


1.167 


2.48 x 10- 


-3 


20.71 


11.82 


9.43 


5.77 


7.71 


1.52 


500.0 


0.5 


0.001 


470.12 


27.60 


2.021 


4.30 x 10- 


-3 


20.71 


11.82 


9.43 


5.77 


7.71 


1.52 



Note - See Note in Table[2]for description. 



power-law. This avoids the need to store the entire star formation 
history and permits more rapid numerical simulations. The time- 
dependent mass return rate from the evolv ing stellar population 
is now computed using the Kroupa] (120051) initial mass function 
and the mass re turn rate as a function of the stellar mass given by 
Maraston (2005), instead of the Ciotti et al. (1991) formulae. For 
the stellar mass-return rate, we also numerically implemented a 
multi-exponential fit. Finally, the dark-matter potential well (and 
the galaxy stellar distribution) are now described as Jaffe (1983) 
models. 

From the astrophysical point of view, in addition to the stan- 
dard outputs considered in SOCS (e.g. final SMBH mass, final 
stellar mass, etc), we now also focus on the duration of the "cold 
phase", of the "obscured phase", and the "Compton thick phase" 
as fiducially computed using the gas temperature and the column 
density of cold gas. The main parameters adopted to fix a model 
are the dark-matter halo mass (distributed to reproduce cosmo- 
logical expectations), the amount of gas deemed to flow onto the 



dark matter halo, and finally the radiative accretion efficiency. 
As a separate parameter, we also consider the initial mass of 
the central SMBH. After some preliminary model exploration 
and verifying that the new treatment, when applied to the SOCS 
galaxy model reproduces the previous results, we performed a 
significant exploration of the parameter space. 

Our main results can be summarized as follows: 

1 . Almost independently of the dark matter halo mass and the 
radiative efficiency, the computed galaxy models have an ini- 
tial phase that extends for some Gyrs, in which the gas tem- 
perature is low and the gas density is high. These galaxies 
would be defined obscured quasars, were we to measure their 
gas column density within an aperture radius of the order of 
fle/100. 

2. At late times, all the SMBHs are found to be in a low ac- 
cretion state without much feedback, and the ISM is overall 
optically thin. The galactic ISM is at about the virial temper- 
ature of the dark matter potential well and should be emit- 



10 



E. Lusso and L. Ciotti: One-zone models for spheroidal galaxies 




Log M„ [M s ] 

Fig. 5. Final SMBH mass versus fi nal stellar mass for all th e ex- 
plored models. The solid line is the Ma rconi & Huntl (l2003) best- 
fit of the Magorrian (1998) relation, and the two dashed lines 
represent the associated l-cr deviation. Different colors indicate 
different initial masses of the M BH , i.e. 10 3 M Q (green), 1O 5 M 
(red, models in Table©, 1O 6 M (blue), and 10 8 M o (black, mod- 
els in Table [2]). Different symbols identify the adopted value for 
radiative efficiency: e = 0.1 (circles), e = 0.01 (triangles), and 
e = 0.001 (squares). The black crosses are the hydrodynamical 
models in Table 1 of Ciotti et al. (2010). 

ting in X-rays, as found for elliptical galaxies in the local 
universe. 

3. Interestingly, we have found that only a specific class of 
models, are found to agree with the observed Magorrian 
relation, i.e. only those with low initial SMBH masses (< 
1O 6 M ) and high radiative accretion efficiency, e ~ 0.1. 
Higher initial SMBH masses, or lower radiative efficiencies 
lead to final SMBH masses that are too high. Therefore, this 
result implies that the seed SMBHs should be quite small, but 
that their mass accretion should occur mainly with high ra- 
diative efficiency, in agreement with observational findings. 

4. In addition, we have also showed how the self-regulated 
Bondi accretion recipe can be easily implemented in numer- 
ical codes, and how it leads to a lower SMBH mass accretion 
than the more common "on-off" Eddington regulation. 

A series of final comments are in order. The first is that the 
presented one-zone models seem to be able to reproduce, with- 
out much fine tuning, two different phases of galaxy evolution, 
namely a first obscured phase where much of the SMBH accre- 
tion and star formation occurs, followed by a hot phase in which 
SMBH accretion is highly sub-Eddington and star formation is 
(almost) entirely absent. The robustness of the two-phase evo- 
lution, characterized both by one-zone and hydrodynamical sim- 
ulations is mainly due to the combined effect of (1) the secular 
decrease in the mass return rate from the evolving stellar popu- 
lations, (2) the time dependence of the SNIa heating (after the 
first Gyrs of evolution), and (3) the cooling function, which are 
identical in the hydrodynamical and the one-zone models. The 



increase in the specific ISM heating with increasing time, and 
that a substantial degassing occurs only when the gas tempera- 
ture is roughly higher by a factor of ~ 4 than the virial tempera- 
ture, leads to the appearance of the two phases in the two types 
of model. 

The second is that to reproduce the present-day Magorrian 
relation, the optimal combination of parameters is a quite low 
initial SMBH mass and a quite high radiative efficiency. We pro- 
pose that this represents a useful constrain of semi-analytical in- 
vestigations (also in the context of high-z galaxy merging). As a 
side-product of this work, we also anticipate that the presented 
multi-exponential fit for stellar evolution and SNIa heating will 
be useful in both hydrodynamical and semi-analytical works. 

As a final comment, we point out a major improvement in 
the code. In a future study, we will relax the assumption of a 
fixed gas density distribution, imposing a time dependence of 
r g (f) as a function of the gas thermal content at that given time. 
We expect that this additional ingredient will cause the models 
to become more sensitive to the adopted parameters on the one 
hand by increasing the gas density during the cold phases (hence 
producing a stronger feedback by decreasing the value of the 
ionization parameter and increasing the mass accretion rate), and 
on the other hand by decreasing the gas density during the hot 
phase, so reducing further the Eddington ratios at late times. We 
naturally, also expect that both the star-formation and accretion 
histories, as well as the durations of the obscured and Compton- 
thick phases to be affected. 
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Appendix A: Interpolating functions for stellar 
evolution 

A. 1 . Stellar mass losses 

The general expression for the stellar mass-return rate from an 
evolving stellar population is 



where 



f 

Jo 



M+(t')W t (t-t')dt', 



(A.l) 



where MX is the instantaneous star-formation rate and W*(t - f') 
is the normalized stellar death rate for a stellar population of age 
t — t' . For an initial mass function (IMF) of unitary total mass 
and age t, 



WJt) = IMF[M T o(f)]|MTo(0|AM w [M T o(f)], 



(A.2) 



where Mxo(0 is the mass of stars ent ering the turn-off at tim e 
t, and AM W their mass loss. Following lCiotfi & Ostrikeri d2007l) . 
we assume that 



( 0.923M-0.48, 0.08 < M < 8.5 
AM w = iM-1.4, 8.5<M<40 
M/2, M > 40. 



(A.3) 



In our models, we adopt a Kroupa (2001) IMF, with a minimum 
mass of 0.08 M Q and a maximum mass of 100 M , while Mjo(t) 
is taken from Maraston (2005) 

log M T0 (f) = 2.982 + 0.213 logf - 0.108(logf) 2 + 
0.006(logr) 3 . 



(A.4) 



In the formula above, t is in yrs and Mjo in solar masses. The 
main computational problem posed by the evaluation of the inte- 
gral in Eq. dA.lt is that, in principle, the entire history of MX(t') 
must be stored, which requires a prohibitively large amount of 
memory and computational time. In the special case of an ex- 
ponential time dependence of W„ this problem can be fortu- 
nately avoided, as shown in ICiotti & Ostrikeri d2001l |2007|) . In 
the present case, W* is not an exponential function, being more 
similar to a power law. However, we assume that a fit in terms of 
exponentials has been obtained, i.e. 







(A.5) 



where the timescales a, and fo, are known. From the above equa- 
tion and Eq. dA.lt . it follows that 

FJf) 



m:(o = v M 



(A.6) 



Fi(t)= f M:(t')e-'^ At' . 
Jo 

We divide the integral into 



(A.7) 



Fi(t) 



J 1-5, 



Mt(t')e At' + 



f 

Jo 



Mt(t')e ~&t', (A.8) 



where 5t represents the last time-step. The evaluation of the first 
integral can be done by using a simple trapezoidal rule and only 
the values of MX at t and t — St are required. The second integral 
is transformed by adding and subtracting 5t in the exponential 
term, and simple algebra finally shows that 



St 



Mt(t) + e ^Mt(t-dt) 



+ e -,Fi{t-5f). (A.9) 



Therefore, the introduction of the multi-exponential fit allow us 
to compute the instantaneous mass return rate at time t just by 
storing the values of the Fj functions at the previous time-step. 
We performed a non-linear fit of the tabulated values of Eq. dA.2t 
(obtained with the exact formulae). We found that an acceptable 
approximation over the entire Hubble time is obtained with the 
coefficients in the first row of Table fA.lK with a percentage error 
of at most of 10% and only at very late times, when the mass 
return rate of a stellar population is negligible). 

A.2. SN la 

In a given simple stellar population, the total (i.e. volume- 
integrated) energy release by SNIa is 



£sn(0 - £sn^sn(0> 



(A. 10) 



where £sn = 10 51 erg is the fiducial energy released by 
SNIa event and /?sn is the SNIa rate. From Eqs. (11)-(12) in 
ICiotti & Ostrikeri d2007), the SNIa rate for a stellar population 
of mass M* can be written as 



«sn(0 



1.6 x 10~ 13 M« / t 



T BQ M \f H 



[yrs '], 



(A.ll) 



where we adopt a stellar mass-to-light ratio Tb© = 5, £>sn = 1, 
h = 0.7, ?h = 13.7 Gyrs, and s = 1.1. In our model with a 
star formation rate of M», it follows that in each time-step an 
amount M*5t of stars is added to the galaxy body. Therefore, the 
instantaneous total rate of SNIa explosion at time t is 



1.6 x 10" 

T B0 M G 



Jo 



t - 1' 



i.i 



df'. 



(A. 12) 



We expanded the dimensionless power-law kernel in the integral 
above again using a (dimensionless) multi-exponential function, 
so that all the considerations in Sect. IA.1I hold and in particu- 
lar Eq. ( IA.9I ). The functions F, now contain new values of b t 
(in units of Gyrs -1 ). The parameters of the multi-exponential fit 
(with percentual errors < 3% over 13.7 Gyrs) are listed in the 
second row of Table TA. II Therefore, 



1.6 x 10 

T B0 M e 



di 



erg s 



(A.13) 



and finally, the average SNIa heating per unit time and unit vol- 
ume, needed in the gas energy equation is 

* 3L S n(0 



iH.SNIa 



87^ 



(A. 14) 



where again we use half-mass values. 
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Table A.l. Multi-exponential expansion parameters (in Gyrs) for stellar mass losses and SNIa rate. Note that for SNIa the coeffi- 
cients a, are dimensionless. 



"(I 


«l 


a 2 


«3 


04 fls 


bo 




b 2 


bi 


b 4 


b 5 


Stellar mass return 


16.38 


1.93 


0.38 






1.93 


0.19 


0.03 








SNIa 


0.18 


3.22 x 1(T 2 


7.15 x 10~ 3 


1.84 x 10~ 3 


5.22 xlO- 4 1.34 x 10~ 4 


7.43 


1.00 


0.22 


5.79 x 10~ 2 


1.72 x 10~ 2 


5.28 x 10~ 3 



A.3. Column density for the Jaffe model 



In our model the gas density at each computation time is as- 
sumed to be proportional to ph and given by eq fl% . The corre- 
sponding projected gas distribution is 



Sg(*) 



JR 



Pg(r)rdr 



(A.15) 



We note that the projected density formally diverges for R — > 0, 
so that, to estimate a realistic value of the column density we 
need to compute a mean value of the projected density within a 
fiducial aperture. The explicit formulae for "L g (R), and the associ- 
ated projected mass encircled with R, Mp g (R), are given in Jaffe 
(1983). In particular, to estimate (Nh) in Eq. (l22l i we consider 



M Pg (R) = 2n f R' ^ g (R')dR' = 

= M g -47T I rp g (r)V^~ 
Jr 

Integration of equation ( IA. 16b yields 



(A. 16) 



R 2 dr 



M Pg (R) 
M„ 



= R 



R arcsechR 



R arcsec/? 



<R < 1; 



(A. 17) 



where R = R/r g and M Pg (r g ) = (tt/2 - l)M g . 



